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ABSTRACT 

The  image  irradiance  equation  constrains  the  relationship  between  surface  orientation 
in  a  scene  and  the  irradiance  of  its  image.  This  equation  requires  detailed  knowledge  of 
both  the  scene  illumination  and  the  reflectance  of  the  surface  material.  For  this  equation 
to  be  used  to  recover  surface  orientation  from  image  irradiance,  additional  constraints  are 
necessary.  The  constraints  usually  employed  require  that  the  recovered  surface  be  smooth. 
We  demonstrate  that  smoothness  is  not  sufficient  for  this  task. 

A  new  formulation  of  shape  from  shading  is  presented  in  which  surface  orientation  is 
related  to  image  irradiance  without  requiring  detailed  knowledge  of  the  scene  illumination,  or 
of  the  albedo  of  the  surface  material.  This  formulation,  which  assumes  isotropic  scattering, 
provides  some  interesting  performance  parallels  to  those  exhibited  by  the  human  visual 
system. 


1  INTRODUCTION 

Most  previous  work  [1-8]  on  the  problem  of  recovering  surface  shape  from  image  shading 
has  been  based  on  solving  the  image  irradiance  equation,  which  relates  the  radiance  of  a 
scene  to  the  irradiance  of  its  image  [1,2]. 1  This  formulation  of  the  relationship  between 
scene  radiance  and  image  irradiance  is  embodied  in  a  first-order  partial  differential  equation 
expressing  scene  depth  as  a  function  of  image  coordinates.  Such  a  formulation  requires 
specific  knowledge  of  not  only  the  reflectance  characteristics  of  the  surfaces  in  the  scene, 
but  also  the  position  and  strength  of  illumination  sources.  The  approaches  to  solving 
this  differential  equation  have  generally  been  either  by  direct  integration  [l]  or  through 
an  iterative  algorithm  that  attempts  to  reduce  the  difference  between  the  predicted  image 
irradiance  and  the  measured  value  [5-7].  Our  interest  is  in  the  iterative  approach  because 
the  alternative  to  it  —  direct  integration  —  requires  specific  boundary  conditions  that  are 
generally  unknown  (in  natural  scenes),  and  its  behavior  when  applied  to  noisy  pictures,  is 
uncertain. 

As  the  image  irradiance  equation  is  a  single  equation  relating  image  irradiance  and  two 
independent  variables  (specifying  surface  orientation),  it  does  not  uniquely  determine  the 
two  independent  variables  for  a  given  value  of  image  irradiance.  Consequently,  when  this 
equation  is  used  to  recover  surface  shape  additional  constraints  are  necessary.  These  may 

1  Image  irradiance  is  the  light  flux  per  unit  area  falling  on  the  image,  i.e.,  incident  flux  density.  Scene 
radiance  is  the  light  flux  per  unit  projected  area  per  unit  solid  angle  emitted  from  the  scene,  i.e.,  emitted 
flux  density  per  unit  solid  angle. 
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be  imposed  by  boundary  conditions,  by  restrictions  on  the  type  of  surface  to  be  recovered, 
or  by  a  combination  of  the  two.  For  some  images,  when  we  can  determine  important 
features  (such  as  the  fact  that  an  edge  is  an  occlusion  boundary  caused  by  a  surface  turning 
smoothly  away  from  the  viewing  direction),  we  can  use  boundary  conditions  to  constrain  the 
solution;  in  large  portions  of  the  image,  however,  we  can  say  something  only  ahout  the  type 
of  surface  we  would  like  to  recover.  To  date  surface  smoothness  is  the  weakest  additional 
assumption  that  has  allowed  surface  shape  to  be  recovered.  Smoothness  normally  signifies 
that  the  surface  is  continuous  and  that  it  is  once  or  twice  differentiahle.  Smoothness,  as  the 
additional  assumption,  has  had  to  play  the  role  of  propagator  of  boundary  conditions  and 
selector  of  the  surface  to  be  recovered.  Is  smoothness  capable  of  these  tasks  in  general  or  is 
its  usefulness  limited  to  special  cases? 

In  the  first  part  of  this  paper  we  describe  the  various  formulations  that  have  employed 
smoothness,  including  a  relaxation  procedure  of  our  own  that  resembles  its  counterpart 
in  engineering;  we  then  present  results  of  our  experiments  with  these  iterative  procedures. 
Assessing  the  usefulness  of  smoothness  in  this  context,  we  conjecture  as  to  its  utility  in  other 
shape- from-shading  formulations. 

Not  all  authors  have  used  smoothness  as  their  additional  constraint;  some  have 
employed  assumptions  about  surface  shape  instead.  The  assumption  that  the  surface  is 
locally  spherical,  i.e.,  that  its  curvature  is  independent  of  direction,  is  strong  enough  to 
allow  but  a  single  interpretation  for  the  surface  orientation,  and  at  the  same  time,  it  is  also 
one  that  enables  recovery  of  the  surface  orientation  by  purely  local  computation  [8].  In 
addition,  this  shape  constraint  eliminates  the  need  to  know  such  parameters  as  illuminant 
direction  and  surface  albedo.2  Assumptions  about  shape  are  being  traded  for  assumptions 
about  reflectance  behavior.  Can  we  formulate  the  shape-from-shading  problem  without- 
having  to  know  the  details  of  the  surface  reflectance  and  without  making  any  assumptions 
about  the  shape  of  the  surface  we  wish  to  recover? 

In  the  new  formulation  presented  in  the  second  part  of  the  paper,  we  assume  that  scene 
materials  scatter  light  isotropically.  We  make  no  assumptions  about  surface  shape  and  we 
do  not  need  to  know  the  parameters  specifying  illuminant  direction,  illuminant  strength, 
and  surface  albedo.  Our  assumptions  are  about  the  properties  of  reflection  in  the  world; 
these  alone  are  sufficient  to  relate  surface  orientation  to  image  irradiance.  In  situations 
in  which  the  assumption  of  isotropic  scattering  is  invalid,  the  formulation  provides  some 
interesting  parallels  to  human  vision. 


2  ITERATIVE  FORMULATIONS  FOR  SURFACE  RECOVERY 

The  image  irradiance  equation  as  presented  by  Horn  [2],  is 

I{x,y)  =  R{p,q)  , 

where  I[x,y )  is  the  image  irradiance  as  a  function  of  the  image  coordinates  *  and  y,  and 
R{p,q)  is  the  surface  radiance  as  a  function  of  p  and  q,  the  derivatives  of  depth  with  respect 


2Surface  albedo  is  the  materia!  reflectance,  i.e.,  the  ratio  of  scene  radiance  to  scene  irradiance. 
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to  the  image  coordinates.  To  derive  this  equation  orthographic  projection  is  assumed;  while 
orthographic  projection  is  inadequate  to  describe  image  formation  it  is  a  good  approximation 
when  the  scene  objects  are  small  compared  to  the  viewing  distance.  In  the  shape-from- 
shading  approach,  it  is  generally  assumed  that  R{p,q)  is  known  for  all  p  and  q  (that  is,  the 
reflectance  map  is  specified).  The  iterative  approach  applies  this  equation  on  a  pixel- by- pixel 
basis,  that  is,  for  pixel  (»,/) 

A'J  =  R{Pi,j ,  Qi,j )  , 

where  /,-j  is  the  image  irradiance  for  (i,/)th  pixel  and  ptJ-,  q^j  is  the  surface  orientation  of 
the  surface  patch  that  is  imaged  at  pixel  For  convenience  we  use  the  notation 

Ri,j  —  R{PiJ,  Qi,j )  • 

If,  at  some  stage  of  the  iterative  procedure,  we  have  assigned  particular  p,-j,  qij  as  the 
surface  orientation  of  the  (i‘,/)th  pixel,  then  the  residual  expression 

ft ,jR  —  ( U,j  ~  Ri,j )2 

specifies  the  error  caused  by  our  assignment  of  surface  orientation.3  If  this  were  our  only 
constraint,  we  could  select  Pij,<?«,y  so  that  %i,jR  —  0-  This  would  guarantee  that  the 
image  irradiance  equation  is  satisfied  pixel  by  pixel,  but,  because  there  are  infinitely  many 
solutions,  we  need  further  constraints  to  reduce  the  number  of  possible  solutions. 

Smoothness  is  usually  introduced  by  specifying  a  relationship  that  we  would  like  to 
have  hold  between  the  surface  orientation  of  the  (i,y)th  pixel  and  its  neighbors.  The  various 
iterative  approaches  [5-7]  differ  in  the  way  this  relationship  is  specified.  Of  course,  at  a 
particular  stage  of  t  he  iterative  process  this  relationship  between  a  pixel  and  its  neighbors 
will  not  be  exact.  Once  again  we  can  specify  a  residual  equation  for  the  error  in  the 
smoothness  relation. 

fi.J  =  [/{.Pi.j  >  Qiji  Pi—  l,j>  ?•—  Qi+lJ  iPi,j—  1 1  9»,  j—  I  >  P  i,j + 1 T  Qi,j+  1  >  ••*)]  > 

where  /  is  the  relationship  between  the  surface  orientation  at  (»,  j)  and  its  neighbors.  An 
example  of  the  type  of  relationship  is  the  difference  between  the  surface  orientation  of  pixel 
(r,  j)  and  the  mean  value  of  the  surface  orientations  of  its  four-neighbors. 

We  have  two  constraints  that  need  to  be  satisfied  simultaneously,  —  one  from  image 
irradiance  and  one  from  surface  smoothness.  At  each  stage  of  the  iterative  process,  the  total 
residual  error  for  pixel  (*,  j)  can  be  described  by 

ti.i  =  KuR  +  b/  , 

where  X  is  a  weighting  factor  that  can  adjust  the  influence  of  the  error  in  image  irradiance 
to  the  error  in  smoothness.4  For  the  image,  the  total  residual  error  is 

c  =  Ee u  ■ 


3The  form  of  the  error  need  not  be  quadratic  —  the  goals  of  such  a  choice  include  simple  final  expressions. 

4 Since  the  error  in  image  irradiance  is  not  necessarily  commensurate  with  that  in  surface  smoothness,  some 
form  of  normalization  is  required. 
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The  allocation  of  surface  orientations  to  all  pixels  should  minimize  this  total  error,  that  is, 


ae 

dpij 


=  0 


ae 

a<7.\y 


=  o 


Differentiating  £  with  respect  to  p,j  and  also  with  respect  to  g,*j  gives  two  equations  for 
each  pixel  in  the  image.  While  complicated  forms  of  the  relationship  between  Pi,j,qi,j  and 
their  neighboring  pixels  will  generally  occur,  we  choose  our  smoothness  relation  so  that  we 
can  arrange  the  equations  in  open  form 


Pi  j  =  Fi{pitj,qi}j,and  p's  and  q's  of  neighboring  pixels)  , 

qij  =  Fzipij,  qij,  and  p's  and  q's  o f  neighboring  pixels)  , 
where  Fi,  and  F%  are  functions. 

We  therefore  have  an  iterative  scheme  that,  given  some  initial  solution,  we  improve  by 
reducing  the  residual  error  in  image  irradiance  and  surface  smoothness.  We  need  to  ask  the 
following  questions  of  such  a  scheme:  Under  what  conditions  will  it  converge  to  a  solution? 
Is  that  solution  unique?  Does  smoothness,  as  defined  by  our  relation,  give  us  the  type  of 
surface  we  want? 


3  SURFACE  ORIENTATION 

There  are  many  equivalent  parameterizations  of  surface  orientation.  Mentioned  pre¬ 
viously  were  the  parameters  p  and  q,  the  derivatives  of  depth  with  respect  to  image  coor¬ 
dinates.  Some  authors  prefer  using  slant  and  tilt  to  specify  surface  orientation.  Slant  is  the 
angle  between  the  surface  normal  and  the  viewing  direction,  while  tilt  is  the  angle  between 
the  image  x  axis  and  the  projection  of  the  surface  normal  onto  the  image  plane.  Other 
parameterizations  [7]  have  been  used  when  particular  properties  of  the  parameterization 
arc  to  be  exploited.  The  parameters  we  use  are  /  and  m:6 

l  —  sin  (7  cos  r  , 
m  =  sin  ersin  r  , 

where  a  is  the  surface  slant  and  t  its  tilt,  l  is  the  component  of  the  surface  normal  in  the 
direction  of  the  x  axis,  and  m  is  the  component  in  the  y  direction.  We  select  this  particular 
parameterization,  as  l  and  m  are  bounded 

0  <  l2  +  m2  <  1  . 

For  surfaces  that  we  can  see, 

0<-<~  , 


6/  and  m  arc  related  to  p  and  q,  l  = 


\/i+p5 


.  . . riL 


+«■* 


'JFFp‘ 
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0  <  r  <  2jt  . 

Consequently  /  and  m  specify  the  surface  norma]  of  an  imaged  surface  without  ambiguity. 


4  FORMULATIONS  USED  FOR  SURFACE  RECOVERY 


To  explore  the  issues  of  convergence,  propagation  of  boundary  conditions,  and  the 
type  of  surface  promoted  by  smoothness,  we  formulate  the  problem  in  two  ways:  one 
that  parallelis  the  technique  previously  described  and,  alternatively,  one  that  resembles  the 
relaxation  method  used  to  solve  structural  engineering  problems. 

The  function  for  scene  radiance,  used  to  create  synthetic  images  for  the  experiment  and 
employed  by  the  shape  recovery  algorithms,  is 

R(l,  m)  =  0.1 569( - - — - - )  +  Max [0.4437 \/l  -  P  -  m2  +  0.3137/  +  0.3137m,  0]. 

This  function  is  appropriate  for  a  scene  that  exhibits  Lambertian  reflectance  and  is  il¬ 
luminated  by  both  a  collimated  source  and  a  uniform  hemispherical  source.  This  illumina¬ 
tion  was  selected  because  it  is  typical  of  the  illumination  of  outdoor  scenes.  The  particular 
numerical  constants  specify  the  light  direction  and  intensity,  and  the  surface  albedo. 

The  first  formulation  is  similar  to  that  described  previously;  we  shall  call  this  the 
‘conventional’  formulation.  From  the  image  irradiance  equation  we  have  the  error  term 


The  smoothness  constraint  is  the  requirement  that  /,j  be  the  average  of  its  four- 
neighbors,  and  that  rriij  be  the  average  of  its  four-neighbors.  The  error  term  for  smoothness 
is 

S  _  n  M-l,j  +  /»+!, j  +  /», j—  l  +  /»,/+!  -.2 

4  ’ 

_l_ ^ ^  TTlj—  i,j  +  m|+i,j-  +  ffljj-i  d~  TTli,j+i  ^2 


Note  that  this  constraint  is  exact  for  a  surface  that  is  spherical. 

Minimizing  £  =  by  differentiating  with  respect  to  /,-y,  and  with 

respect  to  m,j,  and  then  setting  each  result  equal  to  zero,  we  obtain  the  expressions 


m;j  = 


0.4(/i_ij  d"  /i+i,y  d-  /»,y—  i  d-  /i,j+i) 

0.1(/,-ij-i  d-  /.+  i,j+i  +  /,’-i, j  +  i  d-  /»+i,j-i)  — 


0.05(/,_2,j-  d-  /|+2,J  d"  /t,i-2  d"  /t‘J+2)  d- 
dR  I 


0.8\(/,-)j  Ri,j) 


di 


‘.J 


0.4(m,-_ij  d-  mt+l  j  +  mij-i  -I-  m,ii+1)  - 
0.1(rn,_iiJ-_i  d-  mf+jj+i  +  fn,-_iiJ+i  d-  m.+  ij-J  — 
0.05(m,-_2,y  +  ml+2,i  +  mij- 2  d-  m.j+o)  + 


f)  R 

0.8X( 
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We  use  these  (together  with  the  expression  for  as  our  iterative  scheme  to  improve 

on  an  intial  solution. 

The  other  formulation  we  use,  the  ‘engineering’  formulation,  creates  error  terms  from 
the  image  irradiance  equation  and  the  smoothness  constraints,  but  does  not  combine  these 
into  one  term. 

£*\i  =  (Ai  —  Ri,j)  i 

C  S1—U  +  +  thirl  +  N 

4  )  > 

,  5q  ,  +  "»«+  1.J  +  mi,j~ ,  +  ™ij+  ,  x 

fw  =  - 4 - )  • 

We  view  the  £’a  as  residuals  and  apply  the  relaxation  approach,  i.e.,  reduction  of  the  largest 
residuals.  If  £i,jSl  or  £,-jSa  is  selected  for  reduction  we  choose  to  reduce  both,  as  each  is 
independent  of  the  other.  When  is  chosen  for  reduction,  we  do  the  reduction  in  two 

stages  —  one  stage  altering  /,-j  and  the  other  m,-j.  Of  course  we  can  scale  the  residuals, 
reduce  them  from,  say,  the  image  irradiance  equation  to  a  certain  level  before  introducing 
smoothness,  vary  the  amount  of  correction  we  apply,  (e.g.,  we  can  overrelax)  and  the  like. 
In  fact,  we  can  experiment  with  various  relaxation  approaches.  In  this  formulation  major 
changes  in  the  relaxation  scheme  generally  require  minor  programming  changes. 


5  EXPERIMENTAL  RESULTS 

The  test  image  shown  in  Figure  I  is  that  of  a  hemisphere  placed  on  a  plane,  i.e.,  a 
synthetic  image  generated  by  the  reflectance  function  previously  described.  The  collimated 
light  source  is  at  slant  ^  and  tilt  j  —  which  means  that  it  is  at  the  upper  right  as  we 
view  the  image.  We  purposely  avoided  the  case  in  which  the  collimated  source  is  at  the 
same  position  as  the  viewer,  since  the  resulting  symmetric  reflectance  map  might  bias  the 
algorithm  to  return  a  symmetric  surface.  A  synthetic  image  of  a  sphere  was  selected  as 
the  test  image  because  both  the  image  irradiance  equation  and  the  smoothness  relationship 
we  use  hold  exactly.6  The  performance  of  the  algorithm  to  recover  the  surface  shape  could 
be  assessed  without  the  complications  involved  in  using  inexact  models  for  reflectance  and 
smoothness. 

We  need  initial  solutions  to  start  our  iterative/relaxation  procedures.  We  used  four  sets 
of  initial  conditions:  (I)  a  plane  perpendicular  to  the  viewing  direction;  (2)  a  plane  slanted  ~ 
to  the  viewing  direction;  (3)  a  cone  with  its  axis  along  the  viewing  direction;  (4)  the  correct 
solution  perturbed  by  small  random  errors. 

Previous  work  bas  used  boundary  conditions  to  constrain  the  recovered  surface. 
Investigating  this  approach,  we  constrained  the  surface  in  various  ways:  at  the  edge  of  the 
hemisphere,  at  a  closed  curve  lying  on  the  sphere’s  surface,  or  at  individual  points  on  the 
sphere’s  surface.  We  also  used  the  algorithms  without  any  boundary  conditions  whatsoever. 

Since  we  wished  to  investigate  the  extent  to  which  smoothness  could  propagate  bound¬ 
ary  conditions,  we  used  various  image  quantizations,  namely  10  X  16,  32  X  32,  and  64  X  64. 


eThe  smoothness  relationship  does  not  hold  at  the  edge  of  the  hemisphere  where  it  joins  the  plane. 
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The  findings  can  be  characterized  as  follows: 

•  Both  techniques  —  the  engineering  and  the  conventional  method  —  gave  essentially  the 

same  results. 

•  The  engineering  technique  converged  much  faster  than  the  conventional  technique. 

•  Smoothness  propagates  boundary  conditions  by  no  more  than  a  few  pixels 

•  The  initial  solution  largely  predetermines  the  final  one. 

Figures  3-11  display  examples  of  the  results  we  achieved  with  the  conventional  iterative 
scheme;  the  engineering  scheme  gave  essentially  the  same  results.  In  each  of  these  figures 
the  top  left  picture  shows  the  profile  of  the  recovered  surface  (viewed  from  the  bottom 
left  corner)  ,  while  at  the  top  right  we  find  an  image  that  is  the  sine  of  the  surface  slant, 
with  black  representing  0  and  white  1.  The  bottom  left  is  the  cosine  and  the  bottom  right 
the  sine  of  the  surface  tilt,  with  black  representing  -1,  gray  0,  and  white  +1.  The  results 
are  presented  in  this  manner  so  that  the  performance  of  the  algorithms  can  be  evaluated. 
The  profile  can  on  occasion  appear  more  accurate  than  the  individual  surface  orientations 
(as  might  be  expected  of  an  integration  procedure);  at  other  times,  however,  errors  in  the 
surface  orientation  (sometimes  just  from  the  image  quantization)  of  highly  slanted  surfaces 
cause  the  integration  routine  that  produces  the  surface  shape  for  profiling  to  overstate  the 
error.  Figure  2  shows  the  results  that  should  be  obtained  if  the  shape  recovery  algorithms 
recovered  the  surface  exactly. 

Figures  3-6  illustrate  the  effects  of  various  houndary  conditions.  The  errors  at  the  edge 
of  the  sphere  where  it  joins  the  plane  are  expected,  as  smoothness  does  not  hold  there.  Each 
figure  is  the  result  of  320  iterations,  this  being  five  times  the  linear  dimension  of  the  picture 
used.  The  boundary  condition  at  a  point  affects  an  area  of  approximately  10  pixels  in  radius. 
Only  for  Figure  6,  where  a  random  five  percent  of  pixels  were  set  to  their  correct  values,  is 
the  surface  shape  recovered  correctly.  Smoothness  as  a  propagator  affects  but  a  small  area. 
Figures  4,  7,  and  8  illustrate  this  point  further.  Here  various  image  sizes  are  used.  Observe 
that,  as  the  image  size  increases,  the  boundary  conditions  diminish  in  their  effect  and  the 
solution  becomes  progressively  worse.  Figures  4,  9,  and  10  reveal  the  dominant  influence  of 
the  initial  solution.  Figure  11  is  included  to  show  the  effect  of  smoothness  when  X  =  0  — 
namely,  when  image  irradiance  does  not  affect  the  solution  at  all.  This  figure,  obtained  after 
320  iterations,  demonstrates  what  smoothness  alone  can  achieve,  even  when  the  definition 
of  smoothness  is  exact  for  the  viewed  scene  (a  sphere). 

Smoothness  is  a  poor  selector  of  surface  shape  and  a  poor  propagator  of  boundary 
information  when  it  is  used  to  tie  the  surface  orientation  of  a  particular  surface  point  to  those 
of  its  neighbors.  Generally,  in  engineering,  problems  solved  with  relaxation  techniques  are 
formulations  that  relate  a  given  property  at  one  point  to  that  same  property  at  neighboring 
points  by  means  of  differential  relations.  It  is  the  derivative  that  propagates  boundary 
information  and  selects  a  particular  solution  to  be  recovered.  We  present  such  a  formulation 
below  in  an  attempt  to  relieve  smoothness  of  its  role  as  propagator  and  selector. 


6  SURFACE  RADIANCE  AND  ISOTROPIC  SCATTERING 


Our  formulation  of  the  relationship  between  image  irradiance  and  scene  radiance  is 
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I(x,y)  =  R{l,m)  , 


where  I{x,y)  is  the  image  irradiance  at  image  point  x,y  and  R(l,m)  is  the  scene  radiance 
for  a  surface  normal  we  represent  by  l,m.  R  is  a  function  of  the  components  of  the 
surface  normal  and  they,  in  turn,  are  functions  of  image  coordinates.  R(l,m)  specifies 
the  relationship  between  surface  radiance  and  surface  orientation,  while  l{x,y)  and  m{x,y) 
specify  the  relationship  between  surface  orientation  and  image  coordinates.  R(l,  m)  embodies 
knowledge  of  the  nature  of  surface  reflection,  while  l(x,  y)  and  m(x,  y)  embody  the  surface 
shape. 

To  provide  the  additional  constraints  we  need  for  relating  surface  orientation  to  image 
irradiance,  we  introduce  constraints  that  relate  properties  of  R(l,m),  —  that  is,  constraints 
that  specify  the  relationship  between  surface  radiance  and  surface  orientation.  Such  con¬ 
straints  are 

(I  -  l2)Ru  =  (1  -  m2)Rmm  , 

(Rtt  Rmm  )lm  —  (f2  —  TTl  )Rtm  » 

where  Ru  is  the  second  partial  derivative  of  R  with  respect  to  /,  RmTn  is  the  second  partial 
derivative  of  R  with  respect  to  m,  and  Rim  is  the  second  partial  cross-derivative  of  R  with 
respect  to  /  and  m. 

These  two  partial  differential  equations  embody  the  assumption  of  isotropic  scattering 
(Lambertian  reflectance).  For  isotropic  scattering  R{1,  m)  has  the  form 

R(l,  m)  =  al  +  bm  +  c\Jl  —  l2  —  m2  +  d  , 

where  a,b,c,  and  d  are  constants,  their  values  depending  on  illumination  conditions  and 
surface  albedo.  Note  that  l,m,  and  \/l  —  l2  —  m2  are  the  components  of  the  unit  surface 
normal  in  the  directions  x,y,  and  depth.  R(l,m)  can  be  viewed  as  the  dot  product  of 
the  surface  normal  vector  (/,  m,  \/l  —  l2  —  m2)  and  a  vector  (a,  b,c)  denoting  illumination 
conditions.  As  the  value  of  a  dot  product  is  rotationally  independent  of  the  coordinate 
system,  the  scene  radiance  is  independent  of  the  viewing  direction  —  which  is  the  definition 
of  isotropic  scattering. 

It  is  easily  seen  that  R(l,  m)  —  al+  bm  +  c\J  1  —  P  —  m2  +  d  satisfies  the  pair  of  partial 
differential  equations  given  above.  In  the  appendix  we  show  that  R(l,m)  —  al  +  bm  + 
c/T-  W  —  m2  +  d  is  the  solution  of  the  pair  of  partial  differential  equations.  These  partial 
differential  equations  are  an  alternative  definition  of  isotropic  scattering. 

It  is  worthy  of  note  that  R{l,m)  =  al  +  bm  +  c\J  1  —  l2  —  m2  -+•  d  includes  radiance 
functions  for  multiple  and  extended  illumination  sources,  including  that  for  a  hemispherical 
uniform  source  such  as  the  sky.  The  assumption  of  isotropic  scattering  restricts  the  class 
of  material  surfaces  being  considered,  not  the  illumination  conditions. 


7  EQUATIONS  RELATING  SURFACE  ORIENTATION  TO  IMAGE 
IRRADIANCE 

Differentiating 

I{x,y)  =  R(l,m) 

8 


with  respect  to  x  and  y,  we  obtain 


Ix  =  R[lx  +  RmTnx 


Iy  -  ■  R\  ly  +  RmWly  J 

Ixx  ■  Rtl^x  “I"  Rmm^^x  d"  ^Rlm^x^^x  d"  Ri^xx  d~  Rmfflxx  t 

Iyy  R[ily  Rmm^^y  H"  2/?jm /yfTOy  d"  Rl^yy  d~  Rm^^yy  j 

I x y  —z  Rii^x^y  d"  Rmm^x^y  d"  Rlm{^x^y  d"  ty^x  )  d"  Rrfxy  d"  Rm^xy  7 

where  subscripted  variables  denote  partial  differentation  with  respect  to  the  subscript(s). 
From  the  constraints  for  isotropic  scattering,  we  derive  the  relationships 

_  l-  l2 

Rmm  —  Rim 

Substituting  these  relationships  for  Rn  and  Rmm  in  the  expressions  for  Ixx,Iyy,and  Ixy, 
■we  obtain 


['/( 


lm 

1  -  m2 
lm 


)  d"  (  t  )  d-  2 /jfnjjiJ/m  Ixx  Ri^xx  Rm^z 


lm 

2,1  -  f2 


)  d"  ffly  (  )  d"  2 lyffty\Rtm  —  Iyy  Rllyy  Rm^l 


VV 


j  _ j  ^  j2 

[/*/y(  )  d~  TnxWly{  ~j ~  )  d"  Ixttly  d“  ^y^x\Rtm  =  Ixy  Rl^xy  Rm^xy 


lm 


By  removing  Rim  and  substituting  the  expressions  for  Ri  and  Rm,  defined  by  the 
expressions  for  Ix  and  Iy,  we  produce  two  partial  differential  equations  relating  surface 
orientation  to  image  irradiance: 

o6lxx  +  00mxx  a*)lxy  0^ mxy  =  xOIxx  X^i^xy  j 
a6lyy  +  00myy  —  a6lxy  -  06mxy  =  x^yy  ~  xHxy  , 

where 

a  —  Ixmy  —  Iymx  , 

0  ~  lylz  ~  Ixly  t 

7  =  f*2(l  -  m?)+  mx2(l  -  l2)  +  2lxmxlm  , 

6  —  Zy2(l  —  m2)  +  my2(l  —  Z2)  +  2lymylm  , 

6  —  Z*Zy(l  —  m2)  +  mxmy{  1  —  Z2)  +  ( lxmy  +  lymx)lm  , 

X  =  lxmy  -  lymx  . 

These  equations  relate  surface  orientation  to  image  irradiance  by  parameter-free  ex¬ 
pressions.  They  involve  the  derivatives  of  image  irradiance,  but  not  the  image  irradiance 
itself  —  an  important  feature  if  we  conjecture  such  a  model  for  the  human  visual  system. 
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8  RECOVERY  OF  SURFACE  SHAPE  —  A  SPECIAL  CASE: 

A  SPHERICAL  SURFACE 

It  is  difficult  to  solve  the  equations  relating  surface  orientation  to  image  irradiance,  and 
thus  to  recover  surface  shape  from  observed  image  irradiance.  Two  types  of  approaches  are 
possible.  The  two  differential  equations  can  be  integrated  in  a  step-by-step  manner  or,  given 
some  initial  solution,  a  relaxation  procedure  may  be  employed.  The  difficulties  that  arise 
are  two-fold,  numerical  errors  and  multiple  solutions. 

Solutions  of  the  equation  —  0  (the  developable  surfaces,  e.g.,  a  cylinder)  are  also 
solutions  of  the  equations  relating  surface  orientation  to  image  irradiance.  If  the  images 
intensities  were  known  in  analytic  form  then  analytic  solution  of  the  equations  could  employ 
boundary  conditions  to  select  the  appropiate  solution.  However  since  the  analytic  form 
for  the  image  intensities  is  unknown,  numerical  procedures  must  be  employed.  Numerical 
procedures  to  integrate  the  equations  inevidently  introduces  small  errors.  Instability  of  the 
numerical  scheme  seems  responsible  for  such  errors  eventually  dominating  the  recovered 
solution. 

The  alternative,  a  relaxation  procedure  to  solve  the  equations,  has  its  own  difficulties. 
The  difficulties  experienced  in  the  shape-from-shading  methods  discussed  in  the  first  part 
of  this  paper  dictate  caution.  The  importance  of  a  good  initial  solution  for  a  relaxation 
method  cannot  be  overemphasized.  Simplfying  the  two  partial  differential  equations  (using 
additional  assumptions)  provides  a  method  for  obtaining  an  good  initial  solution. 

The  spherical  approximation  assumes  that  we  are  on  a  spherical  surface.  This  implies 
ly  =  0,  mx  =  0,  lx  =  my,  —  namely,  constant  curvature  independent  of  direction.  For  this 
case  the  partial  differential  equations  become  relationships  between  the  second  derivatives 
of  image  irradiance  and  the  components  of  the  surface  normal; 


1  ~  rn2  ./« 

lm  Ixy 


lm  Ixy 

These  results  for  the  spherical  approximation  are  equivalent  to  those  Pentland  was  able 
to  obtain  [8]  through  local  analysis  of  the  surface.  In  addition  to  providing  a  mechanism 
for  obtaining  an  initial  solution  for  a  relaxation-style  algorithm,  their  direct  application 
estimates  the  surface  orientation  by  local  computation  [8]. 

We  are  actively  engaged  in  tbe  development  of  a  relaxation  procedure  to  transform  the 
initial  solution  (given  by  the  spherical  approximation)  into  a  solution  the  satifies  the  full 
equations. 


9  THE  INFLUENCE  OF  BELIEF  ON  THE  PERFORMANCE  OF 
A  VISUAL  SYSTEM 

The  constraints  derived  for  isotropic  scattering  do  not  have  to  be  true  embodiments  of 
the  physical  laws  of  nature;  rather,  they  can  represent  the  beliefs  a  visual  system  possesses 
regarding  those  laws.  In  circumstances  in  which  such  beliefs  do  not  hold,  the  visual  system 
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■will  err  in  predicting  the  world’s  true  nature.  Of  course,  if  the  model  is  not  a  good 
approximation  of  the  physical  laws  of  nature,  the  visual  system  embodying  it  is  useless. 

The  two  constraints  specifying  isotropic  scattering, 

(l-/2)fl,i=(l  —  rn2)Rmm  , 

(Ru  -  Rmm)lm  =  (/2  -  m2)RlTn  , 

obviously  both  hold  when  the  scattering  is  isotropic,  but  what  is  the  situation  for  other 
forms  of  scattering? 

The  images  produced  by  a  scanning  electron  microscope  constitute  an  intriguing  situa¬ 
tion.  The  appropiate  expression  for  scene  radiance  [7]  is 

R(l,m)  =  a(  1+  1  )  , 

\/l  —  P  —  m2 

where  a  is  a  constant.  This  expression  is  quite  unlike  those  for  natural  scenery,  yet  the 
human  visual  system  'sees’  an  image.  Note  that  the  second  constraint  for  the  isotropic 
scattering  case  is  satisfied  by  this  radiance  function,  but  not  the  first.  The  second  constraint 
is  about  surface  tilt,  as  y2l™  2  =  ta^3r,  where  r  is  the  surface  tilt;  the  first  constraint 
introduces  slant.  In  using  the  equations  relating  surface  orientation  to  image  irradiance  to 
recover  surface  orientation,  one  might  expect  them  to  predict  tilt  correctly  for  surfaces  in 
electron  microscope  images,  but  to  err  in  predicting  slant. 

For  other  forms  of  the  scene  radiance  expressions,  neither  constraint  holds.  Specular 
reflectance  has  been  approximated  [2]  by 

R(l,  m)  —  a[6(l  —  l2  —  m2)  +  cl\/ 1  —  l2  —  m2  +  dm\J 1  —  l2  —  m2]n  , 

wrhere  n  is  a  constant,  usually  having  a  value  between  1  and  10  that  determines  the 
‘sharpness’  of  the  specular  peak. 

For  the  maria  of  the  moon,  the  form  of  scene  radiance  normally  used  [2]  is 

a(bl  +  cm 

\J\  —  P  —  m2 

The  constants  a,b,c,  and  d  in  the  above  expressions  are  associated  with  the  strength  and 
position  of  the  light  source,  as  well  as  with  the  surface  albedo. 

The  constraints  do  not  hold  in  either  of  the  preceding  cases.  We  would  expect  a  visual 
system  embodying  them  to  make  errors  under  these  circumstances.  Nevertheless  this  should 
not  induce  us  to  immediately  begin  searching  for  new  visual  beliefs.  After  all  the  human 
visual  system  is  imperfect  under  conditions  of  specular  reflection;  moreover,  people  observed 
the  moon  throughout  history  without  concluding  that  it  was  spherical. 

If  these  constraints  are  incorporated  in  the  human  visual  system,  the  predictions  based 
on  them  —  i.e.,  when  the  visual  system  will  return  ostensibly  ‘correct’  and  ‘incorrect’ 
information  —  could  be  tested  by  psychophysical  experiments.  Such  predictions  together 
with  their  verification  or  refutation  are  being  investigated. 
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10  CONCLUSION 


The  shape-from-shading  task  (recovering  surface  orientation  from  image  irradiance), 
has  meant  finding  a  solution  to  the  image  irradiance  equation.  This  formulation  requires 
that  the  characteristics  of  the  scene  illumination  and  the  surface  material  be  known.  While 
these  requirements  are  difficult  to  satisify, knowing  them  makes  it  possible  to  apply  tbe  image 
irradiance  equation  to  any  scene  material  for  which  the  scene  radiance  function  is  known. 
Such  application,  however,  is  not  without  difficulty,  appropiate  boundary  conditions  are 
needed  and  the  effect  of  image  noise  is  uncertain. 

To  recover  surface  orientation,  relaxation- style  algorithms  based  on  the  image  irradiance 
equation  employ  additional  constraints.  These  constraints,  which  are  needed  to  supplement 
the  uuderdetermined  image  irradiance  equation,  capture  the  concept  of  smoothness.  While 
smoothness  superficially  determines  the  relationship  between  image  irradiance  and  surface 
orientation,  it  is  too  weak  a  concept  to  propagate  boundary  conditions  and  thus  equally 
ineffectual  as  a  means  of  recovering  the  required  solution. 

In  presenting  a  new  formulation  for  the  shape-from-shading  task,  we  have  traded  the 
need  to  know  the  explicit  form  of  the  scene  radiance  function  for  the  assumption  that 
material  scatters  light  isotropically.  This  model  is  applicable  to  natural  scenery  without 
additional  assumptions  about  illumination  conditions  or  the  albedo  of  the  surface  material. 
The  model  also  demonstrates  some  competence  even  when  the  scattering  is  not  isotropic. 
Such  a  model  poses  the  question:  does  the  human  visual  system  embody  a  particular  belief 
about  the  laws  of  scattering  that  it  applies  even  when  these  laws  are  inexact? 

Effective  numerical  procedures  based  on  this  new  formulation  of  the  shape-from-  shad¬ 
ing  task  remain  unknown  and,  are  subjects  for  further  development. 
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APPENDIX 

We  show  that  the  solution  of  the  system  of  partial  differential  equations, 

(1  -  l2)Ru  =  (1  -  m2)Rmm  , 

{Rl l  ~  —  ( / 2  -  m2 )Rtm  , 

where  Ru  is  the  second  partial  derivative  of  R  with  respect  to  /,  Rmm  is  the  second  partial 
derivative  of  R  with  respect  to  m,  and  Rim  is  the  second  partial  cross  derivative  of  R  with 
respect  to  /  and  m,  is 


R(l,  m)  —  al  +  bm  +  c\/l  —  P  —  m2  +  d  , 
where  a,  b,  c,  and  d  are  constants. 

Proof:  Rearranging 

(l-f2)f?„=(l-m2)i?mm  , 

(Ru  -  Rmm)lm  =  (l2  -  m2)Rim  , 


we  obtain 


Rim  = 


Im 


1  -  m2 


Ru  » 


r>  /m  n 

Him  _  ^  -pf'-mm 


Integrating  Rtm  =  jz^Ru  with  respect  to  /,  that  is, 


I  Rim<ll  =  Y^J  tRndl  , 


gives 


Rm  = 


m 


1  —  m2 


(IR,-R)+F1(m)  , 


where  Fx(m)  is  an  arbitary  function  of  m.  Similarly,  integrating  Rim 
respect  to  m  gives 

Rl  =  i  _  p  imRrn  ~  R)+  F?(l)  , 


Rmm  with 
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where  Fn(Z)  is  an  arbitary  function  of  Z.  Rearranging  these  two  equations,  we  get  a  system 
of  two  first  order  partial  differential  equations 

(1  -l2-m2)R,  +  IR  =  ZmF3(m)  +  (1  —  m2)F4(Z)  , 

(1  —  l2  ~  tn2)Rm  +  mR  =  lmF4(l)  +(1  — /2)F3(m)  , 
where  F3(m)  =  (1  —  m2)Fi(m),  and  F4(Z)  =  (1  —  Multiplying  both  equations  by 

the  integrating  factor  (1  —  Z2  —  m2)~a,  we  obtain 

|j[(l  -  Z2  ~  n?)~'R\  =  (1  -  l2  -  m2)-?[/mF3(m)  +  (1  -  m2)F4(Z)]  , 

^[(l-/2-m2)-iF]  =  (l-/2-m2)-?[/mF4(/)  +  (1  -  Z2)F3(m)]  . 

Before  carrying  out  the  integration,  we  can  find  the  form  of  F3(m),  and  F4(Z)  by 
differentiating  the  first  equation  with  respect  to  m  and  the  second  with  respect  to  Z: 

-^-[(1  -Z2- m2)-*/?j  =  (1  -  l2  -  m2)-“[/(l -Z2  +  2m2)F3(m)  +  m(l  -  m2  +  2  Z2)F4(Z) 
alum 

+  /m(l  -  Z2  -  m2)F,3(m)]  , 

-l2-  mz)-*F]  =  (1  —  /2  —  m2)"5  [/( 1  -  Z2  +  2m2 )F3(m)  +  m(l  -  m2  +  2Z2)F4(Z) 

alam 

+  Zm(l -Z2 -m2)F'4(Z)]  , 

where  F'(A:)  indicates  differentiation  with  respect  to  the  independent  variable  k.  Hence, 

F,3(m)  =  F'4(Z)  . 

F3(m)  is  a  function  of  m  and  F4(Z)  is  a  function  of  Z;  this  implies  that 

F3(m)  =  d  , 

F'a  (0  =  <*  , 

where  d  is  a  constant.  Therefore, 

F3(m)  =  dm  +  6  , 

F4(Z)  =  dl  +n  , 

where  a,  and  6  are  constants.  Returning  to  the  integration  step,  we  now  have  the  expressions 
|t[(1  -  Z2  -  m2)-*F]  =  (1  -  Z2  -  m2)"?[Z(6m  +  rf)  +  o(l  -  m2)]  , 

al 

/-[(l-Z2-m2)-  =  F]  =  (l-Z2-m2)-=[rn(a/  +  tZ)  +  6(l-Z2)]  . 
am 

Integrating  the  first  equation  with  respect  to  Z  and  the  second  with  respect  to  m,  we  obtain 
(1  —  Z2  —  —  (6m  +  </)(l  —  Z2  —  m2)~=  +aZ(l  —  Z2  —  m2)-*  +  F6(m)  , 

(1  —  Z2  —  m2)~^R  —  ( aZ  +  <Z)(l  —  Z2  -  m2)~£  +  6m(l  —  l~  —  m2)-^  +  Fe(Z)  , 
where  Fs(m),  and  Fe(Z)  are  arbitary  functions  of  m  and  Z,  respectively.  We  have  two 
expressions  for  R: 

R  =al  +  bm  +  Fs{m)(l  —  l2  —  m2)*  +  d  , 

R  —al  +  bm  +  F3(/)(l  —  Z2  —  m2)*  +  d  , 

which  are  compatible  if 

Fs(m)  =  Fa(Z)  =  c  , 


where  c  is  a  constant.  The  solution  for  R  is 

R  —  al  +  bm  +  c\/l  —  P  —  m2  +  d  . 
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Figure  1  Original  Image. 


Figure  2  'Ideal’  Reeult.  Top  left  -  profile  of  recovered  surface;  top 
right  -  sine  slant,  blacksO,  whitesl;  bottom  left  -  cosine  tilt,  blacks-1, 
gray=0,  white=+l;  bottom  right  -  sine  tilt. 
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Figure  3.  No  boundary  conditions;  planar  initial  solution  perpen¬ 
dicular  to  viewing  direction;  image  quantization  64  X  64. 


Figure  4.  Boundary  at  edge  of  sphere  given;  planar  initial  solution 
perpendicular  to  viewing  direction;  image  quantization  64  X  64. 
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Figure  5.  Boundary  condition  curve  on  sphere’s  surface  (square 
shape);  planar  initial  solution  perpendicular  to  viewing  direction;  image 
quantization  04  X  64. 


Figure  B.  Random  five  percent  of  points  fixed;  planar  initial  solution 
perpendicular  to  viewing  direction;  image  quantization  64  X  64. 
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Figure  7.  Boundary  at  edge  of  sphere  given;  planar  initial  solution 
perpendicular  to  viewing  direction;  image  quantisation  32  X  32. 


Figure  8.  Boundary  at  edge  of  sphere  given;  planar  initial  solution 
perpendicular  to' viewing  direction;  image  quantisation  16  X  16. 
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Figure  0.  Boundary  at  edge  of  Bphere  given;  conical  initial  solution; 
image  quantization  64  X  64. 


Figure  10.  Boundary  at  edge  of  sphere  given;  planar  initial  solution 
slanted  *•  to  viewing  direction;  image  quantization  64  X  64. 
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Figure  11.  Smoothness  constraint  only.  Boundary  at  edge  of  sphere 
given;  planar  initial  solution  perpendicular  to  viewing  direction;  image 
quantization  64  X  64. 
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